#include "volume.hpp"

Volume::Volume() {
  voxels=NULL;
  num=0;
  delta=1.0;
  dim << 0,0,0;
  pos << 0,0,0;
  head << 0,0,0;
  fileName="untitled.zrawx";
  octRoot=NULL;
  enabledVoxelCount=0;
  skeleton=false;
  sparseData=NULL;
}

Volume::Volume(const Vector3i &dim, unsigned char initial) {
  this->dim=dim;
  num=dim.prod();
  delta=1.0;
  pos << 0,0,0;
  head << 0,0,0;
  fileName="untitled.zrawx";
  octRoot=NULL;
  enabledVoxelCount=0;
  skeleton=false;
  sparseData=NULL;
  voxels=new unsigned char[num];
  memset(voxels,initial,num);
}

Volume::Volume(const QString &fileName,bool smallFootprint,bool markInternal) {
  delta=1.0;
  head << 0,0,0;
  octRoot=NULL;
  enabledVoxelCount=0;
  skeleton=false;
  readVolumeFile(fileName);
  sparseData=NULL;
  if (smallFootprint) {
    createSmallFootprint(markInternal);
    delete voxels;
    voxels=NULL;
  } else {
    if (markInternal) markInternalVoxels();
    if (skeleton) findHeadVoxel();
  }
}

Volume::Volume(const Volume &first,const Volume &second) {
  delta=1.0;
  head << 0,0,0;
  octRoot=NULL;
  enabledVoxelCount=0;
  skeleton=false;
  voxels=NULL;
  if (first.size()==second.size()) {
    sparseData=new QHash<long,unsigned char>();
    //first: 2, second: 64
    QSet<long> firstKeys,secondKeys;
    if (first.getSparseData()) {
      firstKeys=first.getSparseData()->keys().toSet();
    } else {
      Volume *firstCopy=new Volume(first);
      firstCopy->createSmallFootprint(false);
      firstKeys=firstCopy->getSparseData()->keys().toSet();
      delete firstCopy;
    }
    if (second.getSparseData()) {
      secondKeys=second.getSparseData()->keys().toSet();
    } else {
      Volume *secondCopy=new Volume(second);
      secondCopy->createSmallFootprint(false);
      secondKeys=secondCopy->getSparseData()->keys().toSet();
      delete secondCopy;
    }
    qDebug() << firstKeys.size() << endl;
    qDebug() << secondKeys.size() << endl;
    QSetIterator<long> i(firstKeys&secondKeys);
    while (i.hasNext()) {
      sparseData->insert(i.next(),1);
    }
    qDebug() << sparseData->size() << endl;
    QSetIterator<long> j(firstKeys-secondKeys);
    while (j.hasNext()) {
      sparseData->insert(j.next(),2);
    }
    qDebug() << sparseData->size() << endl;
    QSetIterator<long> k(secondKeys-firstKeys);
    while (k.hasNext()) {
      sparseData->insert(k.next(),64);
    }
    qDebug() << sparseData->size() << endl;
    dim=first.size();
    num=dim.prod();
    pos=dim.cast<double>()/2;
  } else {
    sparseData=NULL;
    num=0;
    pos << 0,0,0;
  }
}

void Volume::createSmallFootprint(bool markInternal) {
  if (!sparseData) {
    sparseData=new QHash<long,unsigned char>();
    bool lookingForHead=false;
    if (skeleton) lookingForHead=true;
    for (int y=(int)dim.y()-1;y>=0;y--) {
      for (int z=0;z<dim.z();z++) {
        for (int x=0;x<dim.x();x++) {
          long index=z*dim.x()*dim.y()+y*dim.x()+x;
          unsigned char v=voxels[index];
          if (v!=Background) {
            if (markInternal) {
              if (isSurfaceVoxel(x,y,z)) {
                surfaceVoxelCount++;
              } else {
                v|=Internal;
              }
            }
            sparseData->insert(index,v);
            if (lookingForHead) {
              head << x,y,z;
              lookingForHead=false;
            }
          }
        }
      }
    }
    sparseData->squeeze();
  }
}

const QHash<long,unsigned char> * Volume::getSparseData() const {
  return sparseData;
}

Volume::Volume(const Vector3i &dim,unsigned char *data) {
  this->dim=dim;
  delta=1.0;
  num=dim.prod();
  voxels=new unsigned char[num];
  memcpy(voxels,data,num);
  head << -1,-1,-1;
  pos << 0,0,0;
  octRoot=NULL;
  skeleton=false;
  sparseData=NULL;
  enabledVoxelCount=0;
  fileName="untitled.raw";
}

Volume::Volume(const Volume &other) {
  this->dim=other.size();
  delta=other.getVoxelSize();
  num=dim.prod();
  if (other.data()) {
    voxels=new unsigned char[num];
    memcpy(voxels,other.data(),num);
  } else voxels=NULL;
  if (other.getSparseData()) {
    sparseData=new QHash<long,unsigned char>(*(other.getSparseData()));
  } else sparseData=NULL;
  head=other.getHeadVoxel();
  pos=other.position();
  octRoot=NULL;
  skeleton=other.isSkeleton();
  enabledVoxelCount=0;
  fileName="untitled.raw";
}

Volume::Volume(const Slice &slice) {
  dim=slice.size();
  num=dim.prod();
  delta=slice.getDelta();
  pos=-slice.size().cast<double>()*delta/2;
  head << 0,0,0;
  fileName="untitled.zrawx";
  octRoot=NULL;
  enabledVoxelCount=0;
  skeleton=false;
  sparseData=NULL;
  voxels=new unsigned char[num];
  this->slice()=slice;
}

Volume::Volume(const Volume &other,unsigned char *data) {
  this->dim=other.size();
  delta=other.getVoxelSize();
  num=dim.prod();
  voxels=new unsigned char[num];
  memcpy(voxels,data,num);
  sparseData=NULL;
  head << -1,-1,-1;
  pos=other.position();
  octRoot=NULL;
  skeleton=false;
  enabledVoxelCount=0;
  fileName="untitled.raw";
}

Volume::~Volume() {
  if (voxels!=NULL) delete voxels;
  if (sparseData!=NULL) delete sparseData;
  if (octRoot!=NULL) delete octRoot;
}

Slice Volume::slice() {
  return Slice(Vector3i(0,0,0),dim,*this);
}

Slice Volume::slice(const Vector3i &start,const Vector3i &end) {
  return Slice(start,end,*this);
}

const Vector3i & Volume::size() const {
  return dim;
}

const Vector3d & Volume::position() const {
  return pos;
}

void Volume::setPosition(const Vector3d &pos) {
  this->pos=pos;
}

int Volume::size(int index) {
  if (index==0) {
    return dim.x();
  } else if (index==1) {
    return dim.y();
  } else if (index==2) {
    return dim.z();
  }
  return -1;
}

int Volume::maxSize() {
  //unsigned int tmp=dim.x>dim.y?dim.x:dim.y;
  //return tmp>dim.z?tmp:dim.z;
  return dim.maxCoeff();
}

int Volume::minSize() {
  //unsigned int tmp=dim.x<dim.y?dim.x:dim.y;
  //return tmp<dim.z?tmp:dim.z;
  return dim.minCoeff();
}

double Volume::getVoxelSize() const {
  return delta;
}

OctNode * Volume::getOctreeRoot() {
  if (octRoot==NULL) {
    buildOctree();
  }
  return octRoot;
}

void Volume::buildOctree() {
  if (sparseData) return;
  if (dim.x()==dim.y() && dim.x()==dim.z()) {
    int level=(int)log2(dim.x());
    if (1<<level!=dim.x()) return;
    octRoot=new OctNode();
    octRoot->setPosition(0,0,0);
    octRoot->setSize(dim.x());
    octRoot->setLevel(level);
    addSubOctree(octRoot);
  }
}

//TODO check how this works with sparse data
void Volume::addSubOctree(OctNode *subTreeRoot) {
  int level=subTreeRoot->getLevel();
  if (level>1) {
    for (int k=0;k<2;k++) {
      for (int j=0;j<2;j++) {
        for (int i=0;i<2;i++) {
          OctNode *child=new OctNode(subTreeRoot);
          int size=subTreeRoot->getSize()/2;
          Vector3i offset=subTreeRoot->getPosition();
          child->setSize(size);
          //child->setPosition(offset.x+i*size,offset.y+j*size,offset.z+k*size);
          child->setPosition(offset+Vector3i(i,j,k)*size);
          child->setLevel(level-1);
          child->setValue(0);
          subTreeRoot->addChild(child);
          addSubOctree(child);
        }
      }
    }
  } else {
    Vector3i offset=subTreeRoot->getPosition();
    for (int k=0;k<2;k++) {
      for (int j=0;j<2;j++) {
        for (int i=0;i<2;i++) {
          unsigned char value=getVoxel(offset.x()+i,offset.y()+j,offset.z()+k);
          if (value>0) {
            enabledVoxelCount++;
            OctNode *child=new OctNode(subTreeRoot);
            int size=subTreeRoot->getSize()/2;
            child->setSize(size);
            //child->setPosition(offset.x+i*size,offset.y+j*size,offset.z+k*size);
            child->setPosition(offset+Vector3i(i,j,k)*size);
            child->setLevel(level-1);
            child->setValue(value);
            subTreeRoot->addChild(child);
          }
        }
      }
    }
  }
  //prune?
  if (level>1 && subTreeRoot->getNumChildren()>0) {
    QListIterator<OctNode *> i(*(subTreeRoot->getChildren()));
    while (i.hasNext()) {
      OctNode *child=i.next();
      if (child->getNumChildren()==0 && child->getValue()==0) {
        subTreeRoot->removeChild(child);
        delete child;
      }
    }
  }
  if (subTreeRoot->getNumChildren()==8) {
    QListIterator<OctNode *> i(*(subTreeRoot->getChildren()));
    int valueSum=0;
    while (i.hasNext()) {
      OctNode *child=i.next();
      valueSum+=child->getValue();
    }
    if (valueSum==8) {
      subTreeRoot->setValue(1);
      subTreeRoot->removeAllChildren();
    }
  }
}

bool Volume::binaryEqual(const Volume &other) const {
  return sparseData->keys().toSet()==other.getSparseData()->keys().toSet();
}

void Volume::normalise() {
  if (sparseData) {
    QMutableHashIterator<long,unsigned char> i(*sparseData);
    while (i.hasNext()) {
      i.next();
      i.value()=i.value()==0?0:1;
    }
  } else {
    for (int i=0;i<num;i++) voxels[i]=voxels[i]==0?0:1;
  }
}

void Volume::setVoxelSize(double delta) {
  this->delta=delta;
}

int Volume::numVoxels() {
  return num;
}

unsigned char Volume::getVoxel(int index) {
  if (sparseData) {
    if (sparseData->contains(index)) return (*sparseData)[index];
    else return 0;
  } else return voxels[index];
}

unsigned char Volume::getVoxel(int x,int y,int z) {
  long index=z*dim.x()*dim.y()+y*dim.x()+x;
  if (sparseData) {
    if (sparseData->contains(index)) return (*sparseData)[index];
    else return 0;
  } else return voxels[index];
}

unsigned char Volume::getVoxel(const Vector3i &indices) {
  long index=indices.z()*dim.x()*dim.y()+indices.y()*dim.x()+indices.x();
  if (sparseData) {
    if (sparseData->contains(index)) return sparseData->value(index);
    else return 0;
  } else return voxels[index];
}

void Volume::setVoxel(int index,unsigned char value) {
  if (sparseData) {
    if (sparseData->contains(index)) (*sparseData)[index]=value;
    else sparseData->insert(index,value);
  } else voxels[index]=value;
}

void Volume::setVoxel(int x,int y,int z,unsigned char value) {
  long index=z*dim.x()*dim.y()+y*dim.x()+x;
  if (sparseData) {
    if (sparseData->contains(index)) (*sparseData)[index]=value;
    else sparseData->insert(index,value);
  } else voxels[index]=value;
}

void Volume::setVoxel(const Vector3i &indices,unsigned char value) {
  long index=indices.z()*dim.x()*dim.y()+indices.y()*dim.x()+indices.x();
  if (sparseData) {
    if (sparseData->contains(index)) (*sparseData)[index]=value;
    else sparseData->insert(index,value);
  } else voxels[index]=value;
}

bool Volume::isValidIndex(const Vector3i &index) {
  if (index.x()<0) return false;
  if (index.y()<0) return false;
  if (index.z()<0) return false;
  if (index.x()>=dim.x()) return false;
  if (index.y()>=dim.y()) return false;
  if (index.z()>=dim.z()) return false;
  return true;
}

unsigned char * Volume::data() {
  return voxels;
}

const unsigned char * Volume::data() const {
  return voxels;
}

unsigned char * Volume::data(int x,int y, int z) {
  return voxels+z*dim.x()*dim.y()+y*dim.x()+x;
}

const unsigned char * Volume::data(int x,int y, int z) const {
  return voxels+z*dim.x()*dim.y()+y*dim.x()+x;
}

void Volume::setFileName(const QString &fileName) {
  this->fileName=fileName;
}

//void Volume::readVolumeFile(QString fileName) {
//  this->fileName=fileName;
//  QFile file(fileName);
//  QFileInfo fileInfo(file);
//  if (!fileInfo.exists()) {
//    qDebug() << "Error: file does not exist.";
//    exit(1);
//  }
//  file.open(QIODevice::ReadOnly);
//  int fileSize=file.size();
//  if (fileInfo.suffix()=="raw") {
//    QDataStream in(&file);
//    char *buffer=new char[fileSize];
//    dim.x=dim.y=dim.z=(int)ceil(pow(fileSize,1.0/3.0));
//    in.readRawData(buffer,fileSize);
//    voxels=(unsigned char *)buffer;
//    num=fileSize;
//    pos.x=-(dim.x/2.0);
//    pos.y=0.0;
//    pos.z=-(dim.z/2.0);
//  } else if (fileInfo.suffix()=="vox") {
//    QTextStream in(&file);
//    char *buffer=new char[fileSize];
//    in.setIntegerBase(10);
//    int counter=0;
//    char next;
//    while (!in.atEnd()) {
//      in.skipWhiteSpace();
//      in >> next;
//      if (next=='0') {
//        buffer[counter]=0;
//        counter++;
//      } else if (next=='1') {
//        buffer[counter]=1;
//        counter++;
//      }
//    }
//    voxels=new unsigned char[counter];
//    memcpy(voxels,buffer,counter);
//    delete buffer;
//    dim.x=dim.y=dim.z=(int)ceil(pow(counter,1.0/3.0));
//    num=counter;
//    pos.x=-(dim.x/2.0);
//    pos.y=0.0;
//    pos.z=-(dim.z/2.0);
//  } else if (fileInfo.suffix()=="voxml") {
//    VolumeReader reader(fileName);
//    dim=reader.getSize();
//    delta=reader.getDelta();
//    pos=reader.getPosition();
//    num=reader.getNumVoxels();
//    voxels=new unsigned char[num];
//    memcpy(voxels,reader.getVoxels(),num);
//  }
//  file.close();
//}

void Volume::readVolumeFile(const QString &fileName) {
  read(fileName);
}

void Volume::read(const QString &fileName) {
  this->fileName=fileName;
  QFileInfo fileInfo(fileName);
  if (!fileInfo.exists()) {
    qDebug() << "Error: file does not exist.";
    exit(1);
  }
  if (fileInfo.suffix()=="raw") {
    readRawFile(fileName);
  } else if (fileInfo.suffix()=="vox") {
    readVoxFile(fileName);
  } else if (fileInfo.suffix()=="voxml") {
    readVoxmlFile(fileName);
  } else if (fileInfo.suffix()=="pgm") {
    readPgmFile(fileName);
  } else if (fileInfo.suffix()=="rawx") {
    readRawXFile(fileName);
  } else if (fileInfo.suffix()=="zrawx") {
    readRawXFile(fileName,true);
  } else {
    qDebug() << "Error: not a valid format.";
    exit(1);
  }
}

void Volume::writeVolumeFile() { //TODO deprecate
  write(this->fileName);
}

void Volume::write() {
  write(this->fileName);
}

void Volume::writeVolumeFile(const QString &fileName) { //TODO deprecate
  write(fileName);
}

void Volume::recreateVoxels() {
  voxels=new unsigned char[num];
  std::memset(voxels,0,num);
  QHashIterator<long,unsigned char> i(*sparseData);
  while (i.hasNext()) {
    i.next();
    voxels[i.key()]=i.value();
  }
}

void Volume::write(const QDir &dir) {
  qDebug() << this->fileName << endl;
  qDebug() << dir.absoluteFilePath(this->fileName) << endl;
  write(dir.absoluteFilePath(this->fileName));
}

void Volume::write(const QString &fileName) {
  if (voxels==NULL) recreateVoxels();
  QFileInfo fileInfo(fileName);
  if (fileInfo.suffix()=="raw") {
    writeRawFile(fileName);
  } else if (fileInfo.suffix()=="vox") {
    writeVoxFile(fileName);
  } else if (fileInfo.suffix()=="voxml") {
    writeVoxmlFile(fileName);
  } else if (fileInfo.suffix()=="pgm") {
    writePgmFile(fileName);
  } else if (fileInfo.suffix()=="rawx") {
    writeRawXFile(fileName);
  } else if (fileInfo.suffix()=="zrawx") {
    writeRawXFile(fileName,true);
  } else {
    qDebug() << "Error: not a valid format.";
    exit(1);
  }
}

//void Volume::write(QFile &file) {
//  QFileInfo fileInfo(file);
//  file.open(QIODevice::WriteOnly);
//  if (fileInfo.suffix()=="raw") {
//    QDataStream out(&file);
//    char *buffer=(char *)voxels;
//    out.writeRawData(buffer,num);
//  } else if (fileInfo.suffix()=="vox") {
//    QTextStream out(&file);
//    out << QString("%1").arg(voxels[0]);
//    for (int i=1;i<num;i++) {
//      out << QString(" %1").arg(voxels[i]);
//    }
//    out << "\n";
//  } else if (fileInfo.suffix()=="voxml") {
//    QXmlStreamWriter out(&file);
//    out.setAutoFormatting(true);
//    out.setAutoFormattingIndent(2);
//    out.writeStartDocument();
//    out.writeStartElement("voxels");
//    out.writeEmptyElement("size");
//    out.writeAttribute("x",QString::number(dim.x));
//    out.writeAttribute("y",QString::number(dim.y));
//    out.writeAttribute("z",QString::number(dim.z));
//    out.writeAttribute("delta",QString::number(delta,'f'));
//    out.writeEmptyElement("position");
//    out.writeAttribute("x",QString::number(pos.x));
//    out.writeAttribute("y",QString::number(pos.y));
//    out.writeAttribute("z",QString::number(pos.z));
//    out.writeStartElement("data");
//    out.writeAttribute("encoding","base64");
//    QByteArray tmp((char *)voxels,num);
//    QByteArray b64=tmp.toBase64();
//    QString str(b64);
//    int numLines=str.size()/80;
//    int lastLineLength=str.size()%80;
//    int i;
//    for (i=0;i<numLines;i++) {
//      out.writeCharacters("\n    ");
//      QStringRef line(&str,i*80,80);
//      out.writeCharacters(line.toString());
//    }
//    out.writeCharacters("\n    ");
//    QStringRef line(&str,i*80,lastLineLength);
//    out.writeCharacters(line.toString());
//    out.writeCharacters("\n  ");
//    out.writeEndElement(); //data
//    out.writeEndElement(); //voxels
//    out.writeEndDocument();
//  }
//  file.close();
//}

void Volume::readRawFile(const QString &fileName) {
  QFile file(fileName);
  file.open(QIODevice::ReadOnly);
  int fileSize=file.size();
  QDataStream in(&file);
  char *buffer=new char[fileSize];
  dim << 1,1,1;
  dim*=(int)ceil(pow(fileSize,1.0/3.0));
  in.readRawData(buffer,fileSize);
  voxels=(unsigned char *)buffer;
  num=fileSize;
  for (int i=0;i<num;i++) {
    voxels[i]=voxels[i]==0?0:1;
  }
  /*
  pos.x=-(dim.x/2.0);
  pos.y=0.0;
  pos.z=-(dim.z/2.0);
  */
  pos << -(dim.x()/2.0),0.0,-(dim.z()/2.0);
  file.close();
}

void Volume::readVoxFile(const QString &fileName) {
  QFile file(fileName);
  file.open(QIODevice::ReadOnly);
  QTextStream in(&file);
  in.setIntegerBase(10);
  int counter=0;
  char next;
  char *buffer=new char[file.size()];
  while (!in.atEnd()) {
    in.skipWhiteSpace();
    in >> next;
    if (next=='0') {
      buffer[counter]=0;
      counter++;
    } else if (next=='1') {
      buffer[counter]=1;
      counter++;
    }
  }
  voxels=new unsigned char[counter];
  memcpy(voxels,buffer,counter);
  delete buffer;
  dim << 1,1,1;
  dim*=(int)ceil(pow(counter,1.0/3.0));
  num=counter;
  pos << -(dim.x()/2.0),0.0,-(dim.z()/2.0);
  file.close();
}

void Volume::readVoxmlFile(const QString &fileName) {
  VolumeReader reader(fileName);
  dim=reader.getSize();
  delta=reader.getDelta();
  pos=reader.getPosition();
  num=reader.getNumVoxels();
  voxels=new unsigned char[num];
  memcpy(voxels,reader.getVoxels(),num);
}

void Volume::readRawXFile(const QString &fileName,bool compressed) {
  /*
   * steps:
   * 1 - (dimensions) x y z (delta) d (position) x y z
   * 2 - (voxel values) 1 byte per voxel
   */
  QFile file(fileName);
  QFileInfo fileInfo(fileName);
  file.open(QIODevice::ReadOnly);
  QTextStream in(&file);
  //read meta data
  //dimensions
  QString line=in.readLine();
  QStringList list=line.split(QRegExp("\\s+"));
  dim << list[0].toInt(),list[1].toInt(),list[2].toInt();
  //delta
  line=in.readLine();
  delta=line.trimmed().toDouble();
  //position
  line=in.readLine();
  list=line.split(QRegExp("\\s+"));
  pos << list[0].toDouble(),list[1].toDouble(),list[2].toDouble();
  //skeleton
  line=in.readLine();
  skeleton=(line.trimmed().toInt()==1);
  //read voxel data
  num=dim.prod();
  /*
  int filePos=fileInfo.size()-1;
  file.seek(filePos);
  char c;
  file.getChar(&c);
  if (c==10) {
    qDebug() << "File ends with newline.";
    filePos=fileInfo.size()-num;
  } else filePos=fileInfo.size()-num;
  qDebug() << filePos;
  file.seek(filePos);
  */
  file.seek(in.pos());
  QDataStream data(&file);
  char *buffer=new char[num];
  int size=data.readRawData(buffer,num);
  if (compressed) {
    const unsigned char *constBuffer=(unsigned char *)buffer;
    QByteArray unCompressed=qUncompress(constBuffer,size);
    memcpy(buffer,unCompressed.data(),num);
  }
  voxels=(unsigned char *)buffer;
  file.close();
}

void Volume::readPgmFile(const QString &fileName) {
  /*
   * steps:
   * 0 - magic number
   * 1 - dimensions x y z
   * 2 - max value
   * 3 - voxel values
   */
  QFile file(fileName);
  QFileInfo fileInfo(fileName);
  file.open(QIODevice::ReadOnly);
  QTextStream in(&file);
  QString line;
  QRegExp commentRE("^[P#]");
  //read header (skip comments)
  line=in.readLine();
  while (commentRE.indexIn(line)>-1) {
    line=in.readLine();
  }
  //read dimensions
  QStringList list=line.split(QRegExp("\\s+"));
  dim << list[0].toInt(),list[1].toInt(),list[2].toInt();
  pos << -(dim.x()/2.0),0.0,-(dim.z()/2.0);
  //read voxel data
  num=dim.prod();
  int filePos=fileInfo.size()-num;
  file.seek(filePos);
  QDataStream data(&file);
  char *buffer=new char[num];
  data.readRawData(buffer,num);
  voxels=(unsigned char *)buffer;
  for (int i=0;i<num;i++) {
    voxels[i]=voxels[i]==0?0:1;
  }
  file.close();
}

void Volume::writeRawFile(const QString &fileName) {
  QFile file(fileName);
  file.open(QIODevice::WriteOnly);
  QDataStream out(&file);
  unsigned char *buffer=new unsigned char[num];
  memcpy(buffer,voxels,num);
  for (int i=0;i<num;i++) {
    buffer[i]=buffer[i]==0?0:255;
  }
  out.writeRawData((char *)buffer,num);
  delete buffer;
  file.close();
}

void Volume::writeVoxFile(const QString &fileName) {
  QFile file(fileName);
  file.open(QIODevice::WriteOnly);
  QTextStream out(&file);
  out << QString("%1").arg(voxels[0]==0?0:1);
  for (int i=1;i<num;i++) {
    out << QString(" %1").arg(voxels[i]==0?0:1);
  }
  out << "\n";
  file.close();
}

void Volume::writeVoxmlFile(const QString &fileName) {
  QFile file(fileName);
  file.open(QIODevice::WriteOnly);
  QXmlStreamWriter out(&file);
  out.setAutoFormatting(true);
  out.setAutoFormattingIndent(2);
  out.writeStartDocument();
  out.writeStartElement("voxels");
  out.writeEmptyElement("size");
  out.writeAttribute("x",QString::number(dim.x()));
  out.writeAttribute("y",QString::number(dim.y()));
  out.writeAttribute("z",QString::number(dim.z()));
  out.writeAttribute("delta",QString::number(delta,'f'));
  out.writeEmptyElement("position");
  out.writeAttribute("x",QString::number(pos.x()));
  out.writeAttribute("y",QString::number(pos.y()));
  out.writeAttribute("z",QString::number(pos.z()));
  out.writeStartElement("data");
  out.writeAttribute("encoding","base64");
  QByteArray tmp((char *)voxels,num);
  QByteArray b64=tmp.toBase64();
  QString str(b64);
  int numLines=str.size()/80;
  int lastLineLength=str.size()%80;
  int i;
  for (i=0;i<numLines;i++) {
    out.writeCharacters("\n    ");
    QStringRef line(&str,i*80,80);
    out.writeCharacters(line.toString());
  }
  out.writeCharacters("\n    ");
  QStringRef line(&str,i*80,lastLineLength);
  out.writeCharacters(line.toString());
  out.writeCharacters("\n  ");
  out.writeEndElement(); //data
  out.writeEndElement(); //voxels
  out.writeEndDocument();
  file.close();
}

void Volume::writePgmFile(const QString &fileName) {
  QFile file(fileName);
  file.open(QIODevice::WriteOnly);
  QDataStream out(&file);
  QByteArray magicNumber=QString("P5\n").toAscii();
  QByteArray dimensions=QString("%1 %2 %3\n").arg(dim.x()).arg(dim.y()).arg(dim.z()).toAscii();
  QByteArray mmc=QString("255\n").toAscii();
  out.writeRawData(magicNumber.data(),magicNumber.size());
  out.writeRawData(dimensions.data(),dimensions.size());
  out.writeRawData(mmc.data(),mmc.size());
  unsigned char *buffer=new unsigned char[num];
  memcpy(buffer,voxels,num);
  for (int i=0;i<num;i++) {
    buffer[i]=buffer[i]==0?0:255;
  }
  out.writeRawData((char *)buffer,num);
  delete buffer;
  file.close();
}

void Volume::writeRawXFile(const QString &fileName,bool compressed) {
  QFile file(fileName);
  file.open(QIODevice::WriteOnly);
  QDataStream out(&file);
  QByteArray meta=
    QString("%1 %2 %3\n%4\n%5 %6 %7\n%8\n").arg(dim.x()).arg(dim.y()).arg(dim.z()).arg(delta).arg(pos.x()).arg(pos.y()).arg(pos.z()).arg(skeleton?1:0).toAscii();
  out.writeRawData(meta.data(),meta.size());
  /*
  unsigned char *buffer=new unsigned char[num];
  memcpy(buffer,voxels,num);
  for (int i=0;i<num;i++) {
    buffer[i]=buffer[i]==0?0:255;
  }
  */
  if (compressed) {
    QByteArray compressed=qCompress(voxels,num,9);
    out.writeRawData(compressed.data(),compressed.size());
  } else {
    out.writeRawData((char *)voxels,num);
  }
  file.close();
}

int Volume::getNumSurfaceVoxels() {
  return surfaceVoxelCount;
}

int Volume::getNumEnabledVoxels() {
  if (sparseData) return sparseData->size();
  else return enabledVoxelCount; //TODO not reliable!
}

void Volume::markInternalVoxels() {
  enabledVoxelCount=0,surfaceVoxelCount=0;
  for (int z=0;z<dim.z();z++) {
    for (int y=0;y<dim.y();y++) {
      for (int x=0;x<dim.x();x++) {
        if (getVoxel(x,y,z)!=Background) {
          enabledVoxelCount++;
          if (isSurfaceVoxel(x,y,z)) {
            surfaceVoxelCount++;
          } else {
            setVoxel(x,y,z,getVoxel(x,y,z)|Internal);
          }
        }
      }
    }
  }
}

bool Volume::isSurfaceVoxel(int x,int y,int z) {
  if (voxels[z*dim.x()*dim.y()+y*dim.x()+x]) {
    if (x>0) {
      if (voxels[z*dim.x()*dim.y()+y*dim.x()+(x-1)]==Background) return true;
    } else return true;
    if (x<dim.x()-1) {
      if (voxels[z*dim.x()*dim.y()+y*dim.x()+(x+1)]==Background) return true;
    } else return true;
    if (y>0) {
      if (voxels[z*dim.x()*dim.y()+(y-1)*dim.x()+x]==Background) return true;
    } else return true;
    if (y<dim.y()-1) {
      if (voxels[z*dim.x()*dim.y()+(y+1)*dim.x()+x]==Background) return true;
    } else return true;
    if (z>0) {
      if (voxels[(z-1)*dim.x()*dim.y()+y*dim.x()+x]==Background) return true;
    } else return true;
    if (z<dim.z()-1) {
      if (voxels[(z+1)*dim.x()*dim.y()+y*dim.x()+x]==Background) return true;
    } else return true;
  }
  return false;
}

//bogus version that just returns the top-most voxel, replace with face/head detection
void Volume::findHeadVoxel() {
  for (int y=(int)dim.y()-1;y>=0;y--) {
    for (int z=0;z<dim.z();z++) {
      for (int x=0;x<dim.x();x++) {
        if (getVoxel(x,y,z)!=Background) {
          head << x,y,z;
          return;
        }
      }
    }
  }
}

const Vector3i & Volume::getHeadVoxel() const {
  return head;
}

QString Volume::getValidFormats() {
  return QString("Volume files (*.rawx *.zrawx *.raw *.vox *.voxml *.pgm)");
}

QStringList Volume::getValidFormatList() {
  return QStringList() << "*.rawx" << "*.zrawx" << "*.raw" << "*.vox" << "*.voxml" << "*.pgm";
}

bool Volume::isSkeleton() const {
  return skeleton;
}

void Volume::setSkeleton(bool skeleton) {
  this->skeleton=skeleton;
}
